#!/usr/bin/perl
# 2013/11/07 treat new chromosome name
# 2013/07/18 Average enrichment at genes

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if((@ARGV != 4 and @ARGV != 5) or $ARGV[0] eq '--help'){
	die "Usage : $0 -w [wig file] -b [bed file] -c (use category information)\n";
}

my %opt;
getopts("b:w:c:", \%opt);

my $FILE_bed = $opt{b};
my $FILE_wig = $opt{w};


#---------------------------------------
# method
#---------------------------------------
# 1. wig file を読み込む
# 2. bed file を読み込む
# 3. bed file で指定した部位を10%ごとに分割してそれぞれの位置に相当するenrichmentの平均を計算する

#---------------------------------------
# read wig file
#---------------------------------------
my %WIG;
#track type=wiggle_0 name=pKu70
#variableStep chrom=chr1 span=30
#50      33.541859
#80      29.758246
#110     26.788264
#140     24.255730
#170     21.654645
#200     19.378484
{
	my $fh;
	if($FILE_wig =~ /\.gz/){
		$fh = IO::File->new("gzip -dc $FILE_wig |") or die "cannot open $FILE_wig: $!";
	}else{
		$fh = IO::File->new($FILE_wig) or die "cannot open $FILE_wig: $!";
	}
	$fh->getline();
	my $chr = '';
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/chrom=(\w+) /){
			$chr = $1;
		}else{
			my ($pos, $value) = split /\t/;
			if($value < 0){
				$value = 0;
			}
			$WIG{$chr}{$pos} = $value;
		}
	}
	$fh->close();
}


#---------------------------------------
# read bed file
#---------------------------------------
#chr1    1568383 1568732 LTR.31  0       -       1568383 1568732 0       1       349,    0,
#chr1    2173898 2174235 LTR.35  0       -       2173898 2174235 0       1       337,    0,
#chr1    2474234 2474583 LTR.39  0       -       2474234 2474583 0       1       349,    0,
#chr1    2693885 2694225 LTR.41  0       -       2693885 2694225 0       1       340,    0,
#chr1    2843758 2844091 LTR.42  0       -       2843758 2844091 0       1       333,    0,
my %BED;
my %mask;
{
	my $fh = IO::File->new($FILE_bed) or die "cannot open $FILE_bed: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		my ($chr, $start, $end, $name, $group, $direction) = split /\t/;
		$BED{$chr}{$start} = [$end, $direction, $group];
		my $mask_start = &convert($start) - 2000;
		my $mask_end = &convert($end) + 2000;
		for(my $i = $mask_start; $i <= $mask_end; $i += 30){
			$mask{$chr}{$i}++;
		}
	}
	$fh->close();
}

### ２つのBED　で示した位置が近すぎる場合は除く (<2kb以内)
my %BED_far;
{
	foreach my $chr(qw(I II III)){
		foreach my $start(keys %{$BED{$chr}}){
			my ($end, $direction, $group) = @{$BED{$chr}{$start}};
			my $start_c = &convert($start);
			my $end_c = &convert($end);
			my $flag_no = 0;
			for(my $i = $start_c; $i <= $end_c; $i += 30){
				if(exists $mask{$chr}{$i}){
					if($mask{$chr}{$i} > 1){
						$flag_no = 1;
						last;
					}
				}
			}
			if($flag_no == 0){
				$BED_far{$chr}{$start} = [$end, $direction, $group];
			}
		}
	}
}



#---------------------------------------
# calculate enrichment around 2000bp
#---------------------------------------
my %PEAKS;
my %groups;
foreach my $chr(qw(I II III)){
	foreach my $start(keys %{$BED_far{$chr}}){
		my ($end, $direction, $group) = @{$BED_far{$chr}{$start}};
		my $length = $end - $start;
		my $each_length = $length / 10;

		# もし-cオプションを指定していなかったら区別しない
		unless(exists $opt{c}){
			$group = 'average';
		}
		$groups{$group} = 1;

		my $category = 0;
		if($direction eq '+'){
			for(my $i = $start - 2000; $i < $start; $i += 100){
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $i);
				$category++;
			}
			for(my $i = 0; $i < 10; $i++){
				my $location = $start + int($each_length * $i);
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $location);
				$category++;
			}
			for(my $i = $end; $i <= $end + 2000; $i += 100){
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $i);
				$category++;
			}
		}elsif($direction eq '-'){
			for(my $i = $end + 2000; $i > $end; $i -= 100){
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $i);
				$category++;
			}
			for(my $i = 0; $i < 10; $i++){
				my $location = $end -  int($each_length * $i);
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $location);
				$category++;
			}
			for(my $i = $start; $i >= $start - 2000; $i -= 100){
				push @{$PEAKS{$category}{$group}}, &GetWig($chr, $i);
				$category++;
			}
		}
	}
}

my %tick_for_draw = (
	0 => 0,
	1 => 100,
	2 => 200,
	3 => 300,
	4 => 400,
	5 => 500,
	6 => 600,
	7 => 700,
	8 => 800,
	9 => 900,
	10 => 1000,
	11 => 1100,
	12 => 1200,
	13 => 1300,
	14 => 1400,
	15 => 1500,
	16 => 1600,
	17 => 1700,
	18 => 1800,
	19 => 1900,
	20 => 2000,
	21 => 2200,
	22 => 2400,
	23 => 2600,
	24 => 2800,
	25 => 3000,
	26 => 3200,
	27 => 3400,
	28 => 3600,
	29 => 3800,
	30 => 4000,
	31 => 4100,
	32 => 4200,
	33 => 4300,
	34 => 4400,
	35 => 4500,
	36 => 4600,
	37 => 4700,
	38 => 4800,
	39 => 4900,
	40 => 5000,
	41 => 5100,
	42 => 5200,
	43 => 5300,
	44 => 5400,
	45 => 5500,
	46 => 5600,
	47 => 5700,
	48 => 5800,
	49 => 5900,
	50 => 6000
);

my %memori = (
	0 => "-2000bp",
	5 => "-1500bp",
	10 => "-1000bp",
	15 => "-500bp",
	20 => '0',
	22 => '20',
	24 => '40',
	26 => '50',
	28 => '80',
	30 => '100',
	35 => "500bp",
	40 => "1000bp",
	45 => "1500bp",
	50 => "2000bp",
);


my @groups;
foreach my $key(keys %groups){
	push @groups, $key;
}


#---------------------------------------
# output 上下5%ずつのデータは除去して平均を計算する
#---------------------------------------
print "point\ttick\tmemori\t" . join("\t", @groups) . "\n";
foreach my $pos(sort {$a <=> $b} keys %PEAKS){
	my $tick  = $tick_for_draw{$pos};
	my $memori = exists $memori{$pos} ? $memori{$pos} : "";
	print "$pos\t$tick\t$memori";
	foreach my $group(@groups){
		my @sorted = sort {$a <=> $b} @{$PEAKS{$pos}{$group}};
		my $SUM = 0;
		my $NUM = 0;

		my $MAX_index = @sorted - 1;
		my $MIN_index = 0;
		# もしデータ数が10以上あるようなら、上位5%と下位5%を除く
		if(@sorted > 10){
			$MAX_index = int(@sorted * 0.95);
			$MIN_index = int(@sorted * 0.05);
		}

		for(my $i = $MIN_index; $i <= $MAX_index; $i++){
			$SUM += $sorted[$i];
			$NUM++;
		}

		printf "\t%f", $SUM / $NUM;
	}
	print "\n";
}


sub GetWig{
	my ($chr, $pos) = @_;
	for(my $i = 0; $i < 100; $i++){
		my $loc1 = $pos + $i;
		if(exists $WIG{$chr}{$loc1}){
			return $WIG{$chr}{$loc1};
		}
		my $loc2 = $pos - $i;
		if(exists $WIG{$chr}{$loc2}){
			return $WIG{$chr}{$loc2};
		}
	}
	return 0;
}

# マスクするときの計算高速化のため
sub convert{
	my ($num) = @_;
	my $conv = int($num / 30) * 30 + 20;
	return $conv;
}


